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ABSTRACT 

The role played by magnetic fields in the intracluster medium (ICM) of galaxy clusters is complex. 
The weakly collisional nature of the ICM leads to thermal conduction that is channeled along field 
lines. This anisotropic heat conduction profoundly changes the stability of the ICM atmosphere, 
with convective stabilities being driven by temperature gradients of either sign. Here, we employ 
the Athena magnetohydrodynamic code to investigate the local non-linear behavior of the heat-fiux 
driven buoyancy instability (HBI), relevant in the cores of cooling-core clusters where the temperature 
increases with radius. We study a grid of 2-d simulations that span a large range of initial magnetic 
field strengths and numerical resolutions. For very weak initial fields, we recover the previously 
known result that the HBI wraps the field in the horizontal direction thereby shutting off the heat 
fiux. However, we find that simulations which begin with intermediate initial field strengths have a 
qualitatively different behavior, forming HBI-stable filaments that resist field-line wrapping and enable 
sustained vertical conductive heat fiux at a level of 10-25% of the Spitzer value. While astrophysical 
conclusions regarding the role of conduction in cooling cores require detailed global models, our local 
study proves that systems dominated by HBI do not necessarily quench the conductive heat flux. 
Subject headings: conduction - galaxies: clusters: intracluster medium - instabilities - magnetic fields 
- MHD - plasmas 



1. INTRODUCTION 

The intracluster medium (ICM) is a hot, weakly- 
collisional plasma comprising the majority of the bary- 
onic mass in galaxy clusters. It has gained increasing 
attention in the last few decades since it provides a rich 
environment for the study of feedback, cosmology, low- 
density magnetohydrodynamics (MHD), and magnetoge- 
nesis. One area of recent focus is the study of dynamics 
induced by anisotropic conduction and anisotropic vis- 
cosity. The main goal of these studies has been to develop 
a better understanding of the role of magnetic fields in 
moderating heat conduction in the ICM. 

Anisotropic conduction results from the fact that the 
electron gyroradius is much smaller than its mean free 
path and renders the ICM atmosphere buoyantly unsta- 
ble to temperature gradients of either sign, relative to the 
direction of gravity. The outer region of the ICM, where 
virialization supports a negative temperature gradient, 
falls subject to the magnetothermal instability (MTI; 
Balbus 2000). On the other hand, in the inner region 
ot cool-core clusters the plasma density is high enough 
for bremsstrahlung cooling to turn over the temperature 
gradient and the atmosphere is uns table t o the heat-fiux 
driven buoyancy instability (HBI; |Quataert||2008[ here- 
after Q08). 

Previous studies of the local non -linear evolution of 
the HBI using numerical si mulations (Parrish & Quataert' 
20081 iMcCourt et al.|2011|) found that the instability acts 
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to quench vertical heat conduction by wrapping magnetic 
field lines horizontally thereby insulating each tempera- 
ture layer from one another. If this strong insulation ef- 
fect occurs in real clusters, it renders thermal conduction 
irrelevant in consideration of the thermal energy balance 
of the cooling core. 

In the absence of significant conduction, one must in- 
voke AGN feedback or another type of turbulent stirring 
to stave off catastrophic cooling, but any such mecha- 
nism would simultaneously reorient some field lines ver- 
tically and conduction would again drive turbulence via 
the HBI. Therefore any self-consistent model of the ICM 
needs to both realistically capt ure the physics of the HBI 
as well as f eedback processes (Balbus & Reynolds 20^[ 
Sharma et al. 2009; Parrish et al. 2010 , Kunz et al. 20^ 

'I'he HBI is fundamentally a weak-tield instability, op- 
erating even when magnetic forces are completely neg- 
ligible (provided that the field is strong enough to keep 
the electron gyroradius smaller than its mean free path). 
For this reason, most previous studies of the HBI have 
assumed extremely weak initial fields, P > 10^^ where 
P = P/Pmag with P being the thermal plasma pres- 
sure and Pmag the magnetic pressure. However, the role 
of magnetic field/tension on t he non- linear behav ior of 
the HBI is poorly understood. |Kunz et al.| (|2012| here- 
after K2012) performed semi-global HBI runs with ini- 
tial ^ = 10^ and found a qualitatively different behavior 
to that seen previously. Instead of the field line wrap- 
ping and quenching of conduction seen by previous au- 
thors, K2012 found sustained conduction via magnetic 
filaments. However, it was unclear whether the semi- 
global nature of the simulation or the stronger initial 
field was the main cause of this difference. 

This leaves open several questions. Does magnetic ten- 
sion play a defining role in shaping the non-linear behav- 



ior of the HBI? Since tension forces are most relevant 
for highly curved field lines, is the numerical resolution 
typically used in simulations sufficient to fully capture 
the physics of the HBI? What role do magnetic filaments 
play in the transport of heat in the ICM? 

In this paper we present new 2-d simulations in the lo- 
cal regime which cover the resolution-/^ parameter space 
in an attempt to answer these questions. We explore how 
the non-linear dynamics of the weakly collisional ICM 
changes as the initial magnetic field strength is varied 
from complete convective stability (HBI completely sup- 
pressed by magnetic tension at the local scale) to the low 
field strength used in previous studies. In order to de- 
termine whether current simulations have sufficient grid 
resolution for all aspects of the physics of the HBI to 
be fully converged, we also study trends in fundamental 
measures of the non-linear state of the plasma as grid 
resolution is varied. 

We principally find that, as suggested by inviscid sim- 
ulations of K2012, the non-linear state of the HBI is 
qualitatively different than in previous studies when an 
intermediate (within a few orders of magnitude of HBI 
stability) initial magnetic field strength is chosen. In par- 
ticular, there is enough magnetic fiux to form sustained 
vertical filaments which prevent the quenching of vertical 
heat conduction in the ICM, saturat ing at a significant 
fraction of the Spitzer conductivity ( [Spitzer 1962j). We 
present the results of our parameter space exploration 
and describe new insights into the physics of HBI growth 
as well as the formation of magnetic filaments, includ- 
ing a connection between the non-linear turbulent state 
of the plasma and the linear theory. Our results hold 
for moderate resolution 3-d simulations, validating the 
applicability of a 2-d study and the physical model we 
present. 

The paper is organized as follows. In ^ we briefiy 
review the analytic description of the HBIfocusing on 
aspects relevant for our results. In ^ we discuss the 
numerical method, boundary conditions, and initial con- 
ditions, ^presents the results of the simulations and de- 
scribes major trends we observe, which are then placed 
into context of the physics of the HBI in ^ In ^ a 
complimentary 3-d simulation of moderate resolution is 
described, and ^ provides a discussion of our work in 
the context of existing literature. 

2. MHD OF WEAKLY MAGNETIZED PLASMAS 

We use the standard MHD equations with a term 
added to the entropy equation for anisotropic thermal 
conduction along magnetic field lines. For easy compari- 
son to recent literature we primaril y follow the prescrip- 



tion fo r linear analysis laid out in Balbus fc Reynolds] 
(2010, hereafter BRIO) while excluding the term for ra- 
diative cooling. 

We choose rationalized natural units with = c = 
Co = Mo = 1 such that total energy density is given by 
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with the internal (i.e., thermal) energy density e = 
P/(7 — 1), where we assume the adiabatic index 7 = 5/3 
and the mean molecular weight /i = 1 (pure hydrogen, 
fully ionized, plasma). The third term on the right hand 



side is the magnetic energy density, or pressure, used in 
the definition of the plasma /3-parameter 
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The vector quantity x = + defines position in the 
2-d simulations, z points vertically upward and x points 
horizontally to the right. Then, the mass conservation, 
momentum conservation, induction, and entropy equa- 
tions are, respectively. 
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where p is the mass density, v is the fiuid velocity, B is 
the magnetic field vector, P is the gas pressure, g is the 
gravitational acceleration, and Q is the conductive heat 
fiux. D/Dt = d/dt -V is the Lagrangian derivative. 

We define the anisotropic conductive heat fiux Q, with 
b = B/5 the unit vector in the direction of the magnetic 
field, as 

Q = -xb(b ■ V)T, (7) 
where T is the kinetic gas temperature and 
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is the thermal conductivity (Spitzer 1962) in physical 
units. As in Q08 and BRIO we define k = x^/P to be 
the anisotropic thermal diffusion coefficient. In natural 
units, n ^ 10~^ is the Spitzer diffusion coefficient for a 
plasma of temperature 1 keV. 

We ignore the ion component of the heat ffux since 
it is smaller than the electron contributio n by a factor 
of (mi/me)^/^ ^ 42. |Kunz et al.| (|201l|) showed that. 



despite this large ratio, viscosity can have a significant 
effect on the HBI for conductivity at the Spitzer value, 
as chosen here, but that this effect is severely diminished 
deep in cluster cores where the plasma is most collisional. 
For this work, we choose not to include viscosity so that 
we may probe the connection between the formation of 
magnetic filaments and the conduction driven physics of 
the HBI alone. Regardless, the presence of anisotropic 
viscosit y seems to only en hance the formation of fila- 
ments IParrish et aL||2012 ). 



2.1. Background Equilibrium Conditions 

We consider a vertically stratified atmosphere in which 
the gas is not self-gravitating and the extent of the do- 
main is small enough to be in the local regime. We then 
approximate the gravitational acceleration to be a con- 
stant function of position g{z) = —Qqz. 

The initial (unperturbed) magnetic field is purely ver- 
tical for all simulations presented. The simulations start 
with equilibrium initial conditions in which magnetic 
pressure is negligible compared to the thermal pressure. 
Moreover, a uniform initial field provides zero magnetic 
pressure gradient so field strength does not affect the 
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equilibrium condition at all. Therefore the equilibrium 
state of the atmosphere is that of hydrostatic balance 
with 

— = g(^)p = -^op. (8) 

Given an initial temperature profile T(z), in order to 
show distinction between that of the HBI and that of the 
classical adiabatic (Schwarzschild) instability, we employ 
profiles that have entropy increasing with height, i.e., 
that are Schwarzschild stable. Thus, as will be shown 
in ^ we assume a simple power law representation for 
equilibrium temperature and density such that ds/dz > 
where the entropy s is 

_ P 

2.2. Stability analysis 

We refer the reader to BRIO for a full construc- 
tion of the Wentzel-Kramers-Brillouin (WKB) perturba- 
tion analysis using plane wave disturbances of the form 
exp((jt + zk • x) where x is the position vector, a is the 
formal growth rate, and the wavenumber k has Carte- 
sian components. However, there are a few points in the 
analysis which are particularly relevant to the results in 
this paper. 

Without cooling, the linear analysis in BRIO results in 
the following dispersion relation: 

(a + C) (a2+(k-VA)') 
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and the Brunt- Vaisala frequency, N, given by 
, o 1 dP dins 
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describes the buoyant response of an adiabatic plasma. 
k± is the component of k perpendicular to the magnetic 
field direction. The square of the Alfvenic frequency is 

C.i = (k-VA)^ 

This cubic dispersion relation encodes three conditions 
for stabil ity, one of whi ch encodes the HBI and MTI cri- 
terion oflBalbusI (l2000| and Q08, 
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k'^ dz 

where it is convenient to identify the characteristic dy- 
namical timescale with the HBI/MTI growth rate (e- 
folding rate), cj^ = ^ ^l"^ . Eqn. (14) then simplifies. 



in a purely vertical magnetic field 



to 
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This equation is essentially a quantification of the com- 
petition between the effects of a destabilizing tempera- 
ture gradient and the stabilizing influence of magnetic 
tension. 



3. NUMERICAL SETUP 

For our simulatio ns we use the ATHENA MHD code 
(Stone et al. 2008) which uses an unsplit Godunov 
method utilizing constrained transport (Gardiner 



Stone 2005) to preserve the divergence free nature of a 
physical magnetic field. To prevent unphysical transport 
of heat against VT we add a monotonic flux limiter to 
ou r problem alg orithm according to the scheme laid out 
in Sharma & Hammett (2007). 

Atmospheres unstable to the HBI are those with a pos- 
itive temperature gradient. To achieve this condition and 
maintain positive entropy gradient we set up the initial 
equilibrium configuration with temperature, density, and 
pressure given by power laws of the form 
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with Ht = {d\nT/dz) ^ — 2 the characteristic scale 
height of the atmosphere with sound speed Cg = 
^J^P I p ~ 1.3. From hydrostatic equilibrium, with 
^ = 1, we have Tq = 2po = Po = 2/i = 2 and we 
choose an initial magnetic field oriented purely vertically 
(B = 5oi). 

This equilibrium configuration has a conductive heat 
flux that is divergence free. 
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since we choose a constant value for the parallel com- 
ponent of thermal conduction, x\\- Aside from ensuring 
that the HBI will develop with the same growth rate ev- 
erywhere in the anisotropically conducting part of the 
atmosphere, since our models do not account for radia- 
tive cooling and zero-divergence implies no heating, the 
atmosphere is initially in thermal equilibrium. 

Our simulations have dimensionless width x height of 
0.1x0.3 and we limit anisotropic conduction and all our 
measurements of the gas dynamics to the central 0.1x0.1 
region. The remaining 0.1x0.1 buffer zones above and 
below have isotropic conduction of the same magnitude 
as conduction in the active central zone. Early testing 
revealed the vertical reflective momentum boundary con- 
dition to cause suppression of the vertical growth by de- 
flecting motion of the plasma before it naturally reached 
a non-linear state, but addition of buffer zones of this 
size alleviates the effects of the boundaries. 

We choose momentum reflective boundary conditions 
for the upper and lower boundaries to provide mass con- 
servation and pressure support. Periodic boundary con- 
ditions are chosen for the left and right boundaries. The 
temperature is held constant at the upper and lower 
boundaries where the magnetic field vector is set iden- 
tically to that in the last active cell for each column. 
This magnetic boundary condition conserves total verti- 
cal magnetic flux. 

In order to resolve the linear development of the 
HBI cleanly, we applied single-mode perturbations to 
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the equilibrium atmosphere in most of our simulations. 
For consistency, we tested the dependence of our re- 
sults on the perturbation spectrum by running a set 
(the H2dl28_xPS series) of simulations with identical 
initial conditions but a power spectrum perturbation 
in momentum with \vk\^ ex of Kolmogorov type, 

a = 5/3. This produces a non-linear magnetic field 
topology approximately the same as a pure-mode per- 
turbation where two wavelengths fit in each direction in 
the central square box. For all simulations exploring the 
/3-resolution parameter space we start with a single pure- 
mode in both velocity and magnetic field of magnitude 
kx = h= 47r/0.1. 

We performed 24 simulations to cover the /3-resolution 
parameter space and use the designation H2dl28_5, for 
example, to label our 2-d HBI simulation with a resolu- 
tion (covering the active region) of 128 x 128 and an initial 
plasma /3-parameter of /3 = 2 x 10^. To connect our work 
with both linear theory and previous studies we chose five 
values of logio{P) ~ 3,5,7,9 and 11, and five grid res- 
olutions Rx 3R where R = [32,64,128,256,512] In 
addition to these 24 simulations we ran a power spec- 
trum simulation of resolution 128x384 for each P and a 
3-d simulation of resolution 128x128x384 of /3 = 2x10^. 

4. PLASMA BEHAVIOR AS A FUNCTION OF 
RESOLUTION AND /3 

In this section we show the results of our /3-resolution 
parameter space study. We first describe the effect of 
varying /3 (i.e., only changing magnetic field strength) 
on the non-linear saturation of the HBI. We then show 
the dependence of these effects, and evidence for gen- 
eral convergence of the physics of the HBI, with the grid 
resolution. 

As a precursor to our discussion of the full parameter 
space, we show in the first two frames in Figures [l^ and 
the temperature and magnetic field structure of our 
atmosphere under the pure-mode perturbation. These 
frames show the entire height of the atmosphere at times 
when the HBI is in the linear growth phase and just 
reaching the non-linear regime. Note that for both sim- 
ulations of Fig. [T^ and [T]3, corresponding to values of 
f3 = 2 X 10^ and 2 x 10^ respectively, the linear growth 
phase appears identical, in agreement with the expecta- 
tions of the linear theory. 

Now we focus attention on the development of the 
HBI into the nonlinear regime illustrated in the late-time 
frames of Figures [l^ and[T]3. 

4.1. Plasma j3 

Central to any study of the ICM is the plasma {3- 
parameter. In the context of these simulations, consider 
the fully non-linear state of the plasma shown in the later 
frames in Figure [l] The simulation with a moderate ini- 
tial value of /3 ^ 10^ (H2dl28_7) contrasts qualitatively 
to that with p ^ 10^ (H2dl28_9). 

The major qualitative difference is that magnetic fiux 
bundles, or filaments, are present near the end of the 
linear phase in both simulations, but persist over many 

^ We chose not to use limited computing resources on the 25th 
simulation, H2d512_3, since the HBI is stable for that value of /3 
and no dynamics is seen for any grid resolution. 
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Fig. 1. — (a) H2dl28_7 at four times (in units of dynamical times, 
^dy-n)' color bar on the right indicates temperature and mag- 
netic field lines are overlaid, (b) H2dl28_9 at identical times as (a). 
Note the clear distinction between the two after many dynamical 
times, shown in the last frames. 

dynamical times only in the moderate /3 case (2 x 10^ 
and 2 x 10^). These filaments have been seen before as 
the seeds of the horizontal bands of magnetic field lines 
which develop in high p simulat ions. Horizontal mot ion, 
free from any restoring force (see McCourt et al. 2011 for 
discussion), wraps the magnetic field lines thereby insu- 
lating temperature layers from one another and contin- 
uously growing in magnetic field strength. However, for 
the moderate values of j3 these filaments persist against 
horizontal bulk motion and reach an entirely different 
saturated state. 

The most important physical quantity characterizing 
the system, and an excellent discriminator of the different 
non-linear states, is the vertical heat fiux (VHF), 



= - / Xhih- 



V)T • dS, 



where 5* is a surface that cuts horizontally across the do- 
main. Figure [2] shows the horizontal average of the VHF 
across the center of the active domain {z = 1.5) as a func- 
tion of time for all simulations in the pure-mode study. 
The column containing those simulations stable to the 
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Fig. 2. — Horizontally averaged vertical heat flux as a function of time for 24 simulations described in the text. Each panel in the plot 
has identically scaled axes. 



HBI shows constant heat flux consistent with the mag- 
netic field remaining vertical and conducting uniformly. 
The two intermediate (3 columns show a nonzero heat 
flux in the saturated state. As shown in Fig. [l] the HBI 
has shut off VHF in most of the plasma, but the filaments 
remain vertically oriented against the turbulent motion 
of the plasma and support a sustained VHF. The higher 
initial magnetic field strength results in a larger magni- 
tude of sustained VHF, which we explain via the physical 
model for the filaments described in ^ Finally, the high 
13 simulations lead to a non-linear state with zero VHF, 
duplicating the result of previous studies. Thus, one can 
identify the behavior in Fig. 
columns of Figure |2) and Fig. _ 

Another interesting quantity is kinetic energy. In or- 
der to study energy equipartition and turbulence in an 
inherently asymmetric atmosphere we break the kinetic 
energy scalar into the components contributed by mo- 
tion in each direction, KE^ = \'rnv1 and KEx = ^'mvl- 
Figure [s] shows how (KEz), the vertical kinetic energy 
contribution spatially averaged across the entire active 
central region, varies with /3. H2d256_9 and H2d256_ll 
show approximately linear damping of vertical motion, 
but H2d256_5 and H2d256_7 show an increase in vertical 
momentum with higher initial magnetic energy density. 
Both H2d256_5 and H2d256_7 saturate at a value orders 



with the 2nd and 3rd 
with the last two. 
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20 40 60 80 100 120 140 
time (in dynamical times) 

Fig. 3. — The contribution to kinetic energy from vertical motion 
spatiahy averaged across the entire anisotropicahy conducting re- 
gion is shown as a function of time for simulations H2d256_x, x~[5 
(blue-long- dash) , 7 (green- dash- dot- dot- dot) , 9 (black- dash- dot) , 11 
(red- dash)]. 



of magnitude above the weak field cases, where the ver- 
tical momentum does not stop declining up to the end of 
our simulations. 
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However, the kinetic energy is not the only energetic 
discriminant for the non-hnear behavior of the plasma. 
Fig. |4] shows the vertical and horizontal contributions 
to the volume integrated magnetic and kinetic energy 
within the anisotropically conducting region. For j3 = 
2 X 10^ the kinetic energy completely dominates the dy- 
namics of the plasma early on as horizontal motion wraps 
the field into layers. It is only at late times when the 
field wrapping has transferred enough energy into the 
magnetic field for the horizontal component ^^/Stt to 
overcome KE^. 

The moderate f3 columns in Fig. [4] show something 
very different. For sufficient resolution there is approx- 
imate equipartition between BI/Stt ^ KE^ ~ KE^. In 
fact, for initial P = 2 x 10^, even B'^/Stt reaches equiparti- 
tion with the other energy components. In general, how- 
ever, B'^/Stt is subject to a special constraint, namely 
total vertical magnetic ffux through the simulation do- 
main is conserved, and so does not necessarily come into 
equipartition. 

While the simulations start with four different initial 
field strengths spanning four orders of magnitude, the 
total of all the energy components (excluding B'^/Stt) 
is the same for all simulations with resolution at least 
64x192, whether filaments are sustained or not. This is 
expected from the linear theory since the linear growth 
rate of pure-mode perturbations is independent of mag- 
netic field strength (in the weak-field limit) and the linear 
growth entirely determines the amount of kinetic energy 
at the time of transition into the non-linear regime. From 
inspection of Fig. H] it seems that the balance between 
volume integrated kinetic and magnetic energy soon af- 
ter transition into the non-linear regime is an indicator 
of whether or not filaments are sustained. This makes 
sense as an excess in KE is necessary to overcome the 
tension within the filaments. Production of sustained fil- 
aments is intimately tied to the total kinetic energy, and 
thus the turbulent state of the plasma. 

4.2. Grid Resolution 

Now we consider the more subtle effect of grid resolu- 
tion as well as its interplay with p. The amplitude of 
VHF in the two intermediate f3 columns of Fig [2] de- 
creases with increasing resolution. This is illustrated 
in Fig. |5l where the conduction rate vs. time for the 
/3 = 2 X 10^ column at all five grid resolutions is smoothed 
and plotted in superposition. Inset is the time-averaged 
value of conduction during the last 18 dynamical periods 
for each resolution for /3 = 2 x 10^ and 2 x 10^. One can 
see that for /3 ~ 10^ there is negligible change in heat fiux 
for the last doubling of resolution. However, for (3 ^ 10^ 
there is no such sign of convergence. While the data are 
too sparsely sampled in resolution to conclusively make 
this distinction, it is consistent with the general physical 
model of the filaments we propose in the next section. 

Dependence on grid resolution is also evident in Fig. 
|4] Firstly, for moderate /3, increasing resolution produces 
less dispersion in the values of the four volume integrated 
energy components. Secondly, for high /3 simulations 
there is an increase in the vertical energy components 
with increasing resolution. This second effect could be 
due to the resolution of small scale turbulence. 

For high f3 simulations it also takes less time with in- 
creasing resolution for KEx and B'^/Stt to switch dom- 



inance. This is consistent with a lower level of numer- 
ically driven magnetic reconnection in the higher reso- 
lution simulations and the regions of highest magnetic 
field curvature. Finally, in the /3 = 2 x 10^ column, in- 
creasing resolution results in a lower saturated value of 
KE^. In these simulations a significant fraction of the 
vertical motion is bulk fiow along the filament, similar 
to what was seen in the filaments of K2012. We do not 
elaborate on the fiow in this study, but it suffices here to 
simply say the width of filaments is observed to decrease 
with increasing resolution. Thus, the volume integrated 
KEz as well as the resolution dependence of VHF can 
be explained by a dependence of filament width on res- 
olution. In the next section we present a physical model 
for the filaments and explain their dependence on the 
/3-resolution parameter space. 

5. DYNAMICS AND PERSISTENCE OF 
FILAMENTS 

Why do filaments form and what sets their field 
strength and width? We hypothesize that filaments are 
the result of needing to pass a (conserved) vertical mag- 
netic fiux through a plasma subject to an instability 
which attempts to remove vertical field. We suggest that 
the filaments collapse until their internal field strength 
renders them HBI stable, with further collapse likely pre- 
vented by the diffusive process of turbulence. 

Without viscosity, a plasma governed by ideal MHD 
will be driven by turbulent motions toward equipartition 
of the total kinetic and magnetic energies via field am- 
plification. Therefore, after the filaments are produced, 
if the kinetic energy density is sufficiently higher than 
the magnetic energy density, the large scale bulk mo- 
tions of the plasma (preferentially horizontal since the 
atmosphere is classically Schwarzschild stable) will wrap 
up the magnetic field into horizontal layers. If, on the 
other hand, the magnetic field dominates overall, most 
magnetic flux is vertical and the filaments are sustained. 
While the detailed dynamics of filament stability against 
the turbulent motions of the plasma could be framed 
in terms of horizontal velocity shear and height of the 
filaments, the complex magnetic topology and bound- 
ary limitations of the simulations suggest more insight 
may be gleaned from an understanding of filament persis- 
tence in terms of globally averaged quantities. Therefore, 
our model presents a global energetics understanding of 
filament persistence, and a separate model for filament 
width and field strength in terms of the HBI. This separa- 
tion is natural since perturbations of wavelength greater 
than about twice the width of the filament are better un- 
derstood geometrically in terms of bending the fiux tube 
than in terms of the linear HBI. 

A direct prediction of our model is that simulations 
should show that filaments lie on the stable side of the 
threshold between stable and unstable regions of phase 
space. Fig. [6] shows the result of such a measurement, 
obtained via the method we now describe. 

First consider the size of the filaments. Our hypoth- 
esized model predicts filaments will have a width corre- 
sponding to a wavenumber which just satisfies HBI sta- 
bility, i.e., the filaments are stabilized when uj\ > oo'^y^ 
for a region inside the filaments of size equivalent to the 
scale of a random perturbation. 
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Fig. 4. — Evolution of energy components, in our rationalized units, volume integrated over the active (central) region of the computational 
domain of dimensions 0.1x0.1x0.001. The last dimension is the arbitrary skin depth we choose for the compact third dimension. Only 
results of HBI unstable simulations are shown, namely, those for /S = 2 x 10^ are omitted. KEx kinetic energy (black-line), KEz kinetic 
energy {black- dotted) , S^/Stt magnetic energy {red-dash), and ^^/Stt magnetic energy {red- dash- dot). Each panel in the plot has identically 
scaled axes. 
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Fig. 5. — Arbitrarily scaled, horizontallvaveraged vertical heat 
flux plotted in the (3 ~ 10^ column of Fig. [2]versus time in dynam- 
ical units. Colors, in order of transition from black to blue, signify 
increasing resolution. Inset: Each data point represents the flux 
averaged over the last IScjjJ^ for each grid resolution. Conver- 
gence is shown for simulations with initial (3 ~ 10^ (diamonds) but 
not (5 ~ 10'^ (stars). 



In our analysis for Fig. |6) we must be able to define 
and chracterize filaments in an automated manner. We 
consider a filament to be located wherever a^ series of 
horizontally adjacent grid cells satisfy both 16^1 > 0.7 
and Pceii < 0-l(V^)^7 ^ weighted average of P over the 
active domain. The values of magnetic energy (for which 
we can use as a proxy since the average thermal 
pressure varies at most by 2%), the vertical component of 
the temperature gradient, ^\zj and the filament width, 
which we identify with wavenumber , are extracted for 
each such collection of cells for all but the rows of our 
anisotropically conducting active region directly adjacent 
to the buffer regions. Each set of values produces a point 
in Fig. |6] In the plot the value used is calculated 
assuming a wavelength set at twice the measured width of 
the filament. The values of vertical temperature gradient 
within the filament and are unweighted horizontal 
averages across the filament- defining cells. 

If we assume bx is close to zero, which is generally 
observed to be the case for sustained filaments, then it 
follows that K = —k^. If we additionally assume k± ~ 
k\\ for an average perturbation, then 2/cm ~ /c^, and Eqn. 
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Fig. 6. H2d512_5 (upper) and H2d512_7 (lower) filaments. 
Stability of the filaments is demonstrated by plotting the three sta- 
bility diagnostics of Eqn. \19\ measured for each filament. Phase 
space below the solid line is HBI unstable, above is stable. The 
data points are evenly temporally sampled from the last 35.4c(;^J^. 
Those plotted are a subset of all extracted, for visual clarity, but 
the contours represent all the data and enclose 2[^'^' - -'^'^] points. 
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Fig. 7. — Black contours (dash-dot) represent filaments from 
H2dl28_7PS; red (solid) represent those from H2dl28_7. Contour 
levels are set the same as in Fig. [6] 



1 g dT 1 

2T2 /3 



fc^ >0. 



(18) 



This provides a general relation giving the boundary of 
HBI stability in terms of our three diagnostic measures 
of the filaments. In order to facilitate comparison of the 
relation with the data in Fig. [6] we further simplify by 
setting T = 2, in our local simulations where tempera- 
ture does not significantly spatially vary. Also, since the 
width of the filament is the scale at which it is stabilized 
to random perturbations under our model, we identify k 
as the width k..,. We then obtain 



13- 



1 dT^ 

8~dz 



(19) 



In Fig. [6] data points created in the method just de- 
scribed are extracted from the active region for the last 
~ 35 oj'^yn of the simulation at 200 evenly spaced times. 
The magnetic filaments are clearly subject to the bound- 
ary provided by the stability relation, thus supporting 
the theory that the filaments are islands of HBI stability 
within an otherwise HBI susceptible plasma. Specifically, 
consider how the filaments in H2d512_7 and H2d512_5 are 
scattered along the stable side of the curve but generally 
follow it. While there is a significant vertical spread in 
/3~^ above the line, consider that initial values in these 
simulations were = 5 x 10~^ and 5 x 10~^ and the 
final spread in above the intra- filament value is con- 
sistent with the theoretical boundary over nearly three 
orders of magnitude. Note that some filament values of 
I3~^ fall below the initial conditions value. This is consis- 
tent with our definition of a filament having significant 
internal magnetic energy density over the spatial aver- 
age, and most space in the atmosphere has a very weak 
field in the saturated state. 

Figure compares filaments from H2dl28_7 and 
H2dl28_7PS showing that they cover roughly the same 
region in this space. This suggests the pure-mode pertur- 
bation recovers most of the physics of a more multi-mode 
initial condition, and both lie at weaker field strengths 
than the simulations with higher resolution, suggest- 
ing the filaments in the lower resolution simulation are 
prevented from shrinking to obtain their natural peak 
field strength. This picture is also consistent with the 
H2dl28_7 data points residing at the higher values of 
k~'^ on the plot. 

This self-consistent model for the internal stability of 
the filaments is sufficient to explain the width of the fil- 
aments and their dependence on magnetic field strength 
and resolution, as we've shown. However, it also pre- 
dicts the viability of magnetic filaments in a plasma of 
arbitrarily weak magnetic field so long as there is suf- 
ficient resolution for the filament to shrink until it be- 
comes internally HBI stable. This would predict that the 
boundary between our qualitatively different non-linear 
states would shift to higher f3 as resolution is increased in 
Fig. [5] While it is likely we inadequately resolve the /3- 
resolution phase space to see such a trend to begin with, 
the true dissolution of this discrepancy is evident upon 
closer inspection of the global energetics in our simula- 
tions. 

Recall that both the two intermediate and two high P 
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columns of Fig. [2] are well within the weak-field regime 
(See Q. The growth rate then only depends on the 
pure-mode perturbation wavelength, which is identical 
for all simulations in our phase- space suite. Transition 
to non-linearity happens at the same time in all cases, 
and the exponential growth in momentum during the 
linear phase is driven by a conversion of heat into kinetic 
energy. These facts together imply that the total kinetic 
energy imprinted on the plasma at the end of the linear 
regime is the integrated net heat-flux passing into the 
active region during that time. For these reasons, and 
evident in Fig. [4| the peak in kinetic energy at the end 
of linear growth is identical in all the simulations. Those 
with initial /3 ~ 10^ do have a slightly higher transition 
value, but that is consistent with the additional magnetic 
tension maintaining a pure-mode magnetic topology for 
slightly longer, therefore increasing the time over which 
the energy conversion takes place. Clearly the peaks in 
KE components for that column occur at later times than 
in the other columns. 

In the second part of our model we suggest that the 
global persistence of filaments once they've formed de- 
pends on how the imprinted kinetic energy just described 
compares to the initial total magnetic energy. We argue 
that if approximately twice the magnetic energy density 
is greater than the kinetic energy density imprinted at 
the transition to non-linearity, then sustained vertical 
filaments of the form we see will be present. This is ef- 
fectively a statement that either the magnetic field or 
the turbulent motion dominates and one or the other 
wins out early in our simulations, leading to a distinct 
division in parameter space. The 'twice' comes from al- 
lowing an x-component of the filaments to grow to equal 
magnitude with the vertical z-component. Beyond that, 
the filament would be more horizontal than vertical, sig- 
nificantly quench vertical heat flux, and effectively be 
stretched out beyond our definition of a filament. 

This scale-invariant method for understanding filament 
stability is particularly suited to a local study such as 
ours. Of course, it suggests that a sustained source of 
turbulence may significantly alter the non-linear state 
of the plasma. We leave a discussion of this and other 
implications for the last section. 

6. 3-D SIMULATION 



Some studies suggest (Parris h et al.|2Q12 ) that the fila- 
ments formed in 2-d simulations are a result of field lines 
not being able to slip past one another, that is, a pre- 
vention of flux- interchange modes in the plasma. K2012 
comments on this in their H3dBrag simulation, where 
the filaments become more reoriented (horizontally) in 
the lower part of the atmosphere where the effects of 
Braginskii viscosity are less pronounced. However, this 
region of HSdBrag is very close to the simulation bound- 
ary, and the effect of the boundary on the detailed field 
structure is unclear. 

Our model explaining the presence of the filaments, 
presented in the previous section, suggests this existence 
is rooted in the physics of the HBI and not a feature of 
dimensionality. K2012 does not provide a complimentary 
3-d semi-global simulation of the HBI without viscosity 
and so in extension of that work, and to confirm the 
robustness of our results, we performed a 3-d simulation 
of 128x128x384 resolution and moderate field strength. 



streamline 
Var: Speed 
^_-Q.Q4551 

-0.03414 

m- 0.02276 

-0,01139 
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Max. 0,04551 
Min: l,334e-05 




Fig. 8. — Shown is a frame from the saturated state of 
H3dl28_7PS. Color value indicates magnetic field strength and field 
lines are integrated from grids of points at the top and bottom of 
the domain. Many of the field lines thread a filament near the 
center of the anisotropically conducting region. 



H3dl28_7PS. The simulation was in every way identical 
to H2dl28_7PS save for the added dimension. 

Figure [8] shows the magnetic structure of the non-linear 
state of H3dl28_7PS at late times, with coherent fila- 
ments of a qualitatively similar sort to those found in 
the 2-d analog. Just as is pointed out in K2012, these 
filaments contain more magnetic flux because it is eas- 
ier to collect with the additional dimension during linear 
growth. 

The VHF measured in the same way as before is shown 
in Fig. [9] H3dl28_7PS clearly runs long enough for the 
VHF to stabilize at a value of ^ 25% the initial value, just 
as before. This value is no coincidence. The addition of a 
third dimension does not change the internally stabilizing 
P of the filaments, and so the ratio of conducting surface 
area to non-conducting is still set by the initial plasma /3 
and total vertical flux conservation, assuming that most 
vertical flux is bundled into the filaments. 

7. DISCUSSION AND CONCLUSIONS 

In this paper, we present a systematic investigation of 
the local non-linear behavior of the HBI as a function 
of initial magnetic field strength and numerical resolu- 
tion. We identify a previously unrecognized regime of 
behavior at intermediate field strengths in which mag- 
netic filaments are formed that can resist the horizontal 
field line wrapp ing seen in weak- field simulations (Par-] 
Irish fc Quataert 2008; McCourt et al.| |2011[ ) and lead to 
sustained conductive heat flux. We show that the fila- 
ments themselves are HBI-stable and provide the main 
conduit by which (conserved) vertical magnetic flux is 
passed through the domain. In addition to this model for 
internal stabilization, we describe how global energetics 
determine the ultimate qualitative state of the non-linear 
regime, that is, whether or not filaments persist and there 
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Fig. 9. — Vertical Heat Flux integrated over the horizontal plane 
slicing through the center of the active domain in H3dl28_7PS. 



is a resulting significant vertical heat flux. 

The robust formation of filamentary structures is tan- 
talizing since we see this type of structure in the cooling 
cores of galaxy clusters, namely in the form of Ha fil- 
aments. For instance, HST observations of the Perseus 
cluster, which has been studied closely for some time 
( Lynds| p^97Q), reveal complex braided filaments which 
can be as narrow as 50 pc but may be kiloparsecs in 



5( 



length ( Fabia n et al. 2QQ8|). M ore recent observations 



by [McDonald et al. 1 [2010[ 2011) have shown that essen 



tially all cool-core clusters exhibit extended Ra systems 
which appear to extend out to almost, but never beyond, 
the cooling radius of the X-ray emitting ICM. This con- 
nection between the cooling radius and the filamentary 
structure suggests they may condense out of the ICM as 
a result of thermal instability (Field 1965). 

However, it is the local thermal instability which pro- 
duces cold filaments in simulations, with anisotropic con- 
duction primarily responsible for their filamentary struc- 
ture. Without a heating source to offset the background 
global thermal instability, namely, when net cooling is 
purely a function of local density and temperature of the 
gas, buoyancy effec ts stabilize the plasm a to the local 
thermal instability ( |Balbus fc Soker|1989D . K2012 shows 
that for an atmosphere with intermediate magnetic field 
strength and anisotropic viscosity, viscous heating dur- 
ing formation of the filaments and a sustained vertical 
heat flux can significantly delay the cooling catastrophe 
resulting from global thermal instability. This allows lo- 
cal thermal instability to develop and cold filaments to 
naturally form. Our study supports K2012, firmly es- 
tablishing the robust formation of filaments over a range 
in magnetic field strength relevant to cool-core clusters, 
and exploring in detail the magnitude of vertical heat 
flux sustained by the atmosphere and its dependence on 
local magnetic field strength and numerical resolution. 

Though our work has not yet shown a definitive con- 
nection to the observed Ha filaments, in part because 
we do not include radiative cooling in our simulations. 



our findings strongly motivate further work in develop- 
ing physically realistic global simulations of cluster cores, 
where the ICM is likely HBI unstable. Since filaments 
provide a mechanism for transporting heat (and possibly 
cold gas) radially in these atmospheres, an understand- 
ing of the magnetic topology from the local to the global 
scale and how these scales connect will be instrumen- 
tal in understanding the role of conduction in regulating 
global thermal stability. Also instrumental is the inter- 
play between the filaments and astrophysical turbulence, 
since turbulence can ultimately destroy otherwise stable 
filaments. Our discussion of the presence of filaments 
as a result of total energy balance provides some insight 
into the level of turbulence allowed when filaments are 
present, but more testing, especially on the semi-global 
and global scales is required. 

We find, however, that to fully resolve the physics of 
the non-linear HBI in the local regime requires a high 
grid resolution. Thus, feasible global simulations of the 
ICM may require semi-analytic prescriptions for the dy- 
namics of anisotropic conduction as well as anisotropic 
viscosity. This study may aid in production of those 
models since filaments have a key role in the thermody- 
namics of anisotropically conducting plasma. 

Another motivation for this work was to find, if any 
exist, robust measures for convergence with resolution. 
Given the intimate connection between the topology of 
the filaments and more global quantities such as heat con- 
duction, it is not surprising that we find that adequate 
resolution of the dynamics of the filaments is needed in 
order to reach convergence in the plasma diagnostics such 
as VHF and equipartition. The most stringent measure 
we find is to measure the position of the filaments in HBI 
stability space presented in Fig. [6) However, a much sim- 
pler measure that also performs well is to see how much 
the VHF changes with increasing resolution. 

Should the HBI be driving dynamics in relatively large 
regions of real clusters, numerical studies including this 
one suggest that it may have a key role in cluster thermo- 
dynamics. While there is clearly an interaction between 
AGN and the ICM (e.g. radio lobes, etc.), it is not known 
exactly how this feedback may couple effectively to the 
thermodynamics. This work and K2012 suggest there 
may still be a potential for conductive solutions to the 
cooling flow problem. Certainly at this point, the field is 
still in the exploration stage of the fundamental physical 
processes which work together to govern the ICM. 
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